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The nonlinear dynamical behavior of a micromechanical resonator acting as one of the mirrors in 
an optical resonance cavity is investigated. The mechanical motion is coupled to the optical power 
circulating inside the cavity both directly through the radiation pressure and indirectly through 
heating that gives rise to a frequency shift in the mechanical resonance and to thermal deformation. 
The energy stored in the optical cavity is assumed to follow the mirror displacement without any 
lag. In contrast, a finite thermal relaxation rate introduces retardation effects into the mechanical 
equation of motion through temperature dependent terms. Using a combined harmonic balance and 
averaging technique, slow envelope evolution equations are derived. In the limit of small mechanical 
vibrations, the micromechanical system can be described as a nonlinear DufBng-like oscillator. Cou- 
pling to the optical cavity is shown to introduce corrections to the linear dissipation, the nonlinear 
dissipation and the nonlinear elastic constants of the micromechanical mirror. The magnitude and 
the sign of these corrections depend on the exact position of the mirror and on the optical power 
incident on the cavity. In particular, the effective linear dissipation can become negative, causing 
self-excited mechanical oscillations to occur as a result of either a subcritical or supercritical Hopf 
bifurcation. The full slow envelope evolution equations are used to derive the amplitudes and the 
corresponding oscillation frequencies of different limit cycles, and the bifurcation behavior is ana- 
lyzed in detail. Finally, the theoretical results are compared to numerical simulations using realistic 
values of various physical parameters, showing a very good correspondence. 



I. INTRODUCTION 

The experimental study of interactions between light 
and mechanical systems was pioneered more than a hun- 
dred years ago by Crookes Lebedew Q and others 
The two main coupling mechanisms between radia- 
tion and mechanical systems, namely, radiation pressure 
and thermal effects, were already present in these first 
experiments. Since then, the effects of radiation pres- 
sure have attracted a significant interest. An early ex- 
ample is the proposition to use the radiation pressure 
as a driving force in space Another example comes 
from the efforts to detect gravitational waves. The op- 
tomechanical coupling as a source of additional noise in 
gravitational waves detectors and the possibility to utilize 
a high-finesse optomechanical cavity for noise reduction 
in these detectors has been actively discussed for several 
decades (see Refs. [5-9] and references therein). More re- 
cently, similar mechanical mode cooling techniques based 
on radiation pressure have been proposed as a possible 
way to quench the thermal noise in a single mechanical 
vibration mode down to the quantum limit (9|-[T^. 

The renormalization of the effective mechanical damp- 
ing due to coupling of a mechanical oscillator to an op- 
tical resonance cavity is at the heart of these "cooling" 
methods. The root cause of the changes in the effec- 
tive mechanical dissipation in optomechanical systems is 
the retardation in the radiation induced forces. In many 
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studies, a retardation which occurs in the radiationpres- 
sure in optomechanical cavities with high finesse 0, [isl - 
fisl is considered. In such cavities, the optical relaxation 
rate is comparable to the period of the mechanical os- 
cillations. However, high finesse cavities require state of 
the art manufacturing technology and, in general, are 
not readily adjustable for a wide range of different me- 
chanical mirrors. On the other hand, optically induced 
thermal effects have been shown experimentally to affect 
the dynamics of optomechanical systems, including those 
with finesse of order of unity [18!J23|. In these cases, the 
retardation is due to a finite thermal relaxation rate 124- 




In contrast with the thoroughly investigated mechani- 
cal mode cooling effect, other dynamical phenomena that 
arise from the optomechanical coupling have received lim- 
ited theoretical attention. These phenomena include self- 
excited oscillations (l^. l20l [23. 24. 27-31], and changes 
in the effective nonlinear elastic and dissipative behavior 
of an optomechanical mirror. 

As the field of nano optoelectromechanical systems 
(NOEMS) (3^ - [3i ] grows and matures, and, in parallel, 
the search for mechanical systems at quantum limit in- 
tensifies, an increasing number of different optomechan- 
ical systems are being investigated. A theoretical model 
that accurately describes all the phenomena in an op- 
tomechanical system and which is able to reproduce the 
transitional dynamics as well as the steady state and the 
small vibrations behavior would be of great benefit, espe- 
cially for the design of such systems and the experimental 
identification of their parameters. 

In this work, we develop a theoretical model of a mi- 



2 



cromechanical mirror acting as a part of an optical res- 
onance cavity. The mirror is described as a nonlinear 
oscillator, with cubic elastic and dissipative terms in its 
equation of motion f29|, [H, [s^ . The forces acting on the 
mirror include direct radiation pressure, a thermal force 
proportional to the temperature change of the mirror, 
and an external excitation. In addition, a linear depen- 
dence of the mechanical resonance frequency on the tem- 
perature is assumed. U sing a combined harmonic balance 
and averaging method |37| to solve the weakly nonlinear 
equations of motion, we find a practical approximation 
of this model in the form of evolution equations that de- 
scribe the slow envelope dynamics of the system. We 
investigate two important limiting cases of these general 
evolution equations. 

First, we derive the evolution equations for the case of 
small vibrations. In addition to the renormalization of 
the linear mechanical dissipation, we find that the cou- 
pling to an optical resonance cavity introduces additional 
elastic and dissipative nonlinearities into the dynamics 
of the micromechanical mirror. Based on these results, 
stability criteria are derived for small oscillations of the 
mirror, and are shown to coincide with the predictions of 
a local stability analysis of the full dynamical system. In 
addition, the small limit cycle amplitude and frequency is 
given for cases in which a supercritical Hopf bifurcation 
occurs, and the divergence time scale is estimated for 
a stability loss process that leads to a subcritical Hopf 
bifurcation and, consequently, to a jump to a large am- 
plitude limit cycle. 

Next, we explore the behavior of the system at finite 
amplitudes without external excitation. Using the full 
slow envelope evolution equations, we derive the expres- 
sions governing the amplitudes and frequencies of all limit 
cycles fss'l that exist in the system. The resulting steady 
state amplitude equations have the same form as those 
derived in literature from general power or force balance 
considerations (23i. [29l. [sol. |40| . However, in this work, we 
are able to formulate the full evolution equations. There- 
fore, the dynamics of the system can be traced in time, 
in addition to the final steady state solutions similar to 
those previously given in the literature. 

Finally, we explore the validity of our combined har- 
monic balance - averaging method and other assump- 
tions. We find that the method is applicable to a 
wide range of practical optomechanical cavities, espe- 
cially those in which the finesse is relatively low, the 
mechanical quality factor is large and the dependence 
of the mechanical frequency on radiation heating is rel- 
atively weak. In contrast, the amplitude of the mechan- 
ical mirror vibration does not have to be small, and 
can be comparable to the optical wavelength or larger. 
These assumptions are correct for most optomechanical 
resonators, except for those designed specifically to be 
incorporated in high finesse optical cavities. However, 
the mathematical method described here can be readily 
applied to these systems as well. 

In order to experimentally validate the theoretical re- 



sults derived in this article, we have recently studied an 
optomechanical cavity with a moving mirror in the form 
of a freely suspended micromechanical resonator. Using 
the theoretical model developed here, we have been able 
to quantitatively describe the dynamics of micromechan- 
ical mirrors with two different geometries and material 
compositions (4l| . The theory and the experiment have 
been found to be in a good agreement both in the domain 
of forced oscillations and self excitation. 



II. THEORETICAL MODEL 
A. Optomechanical resonance cavity 

Consider an optical resonance cavity constantly 
pumped by monochromatic laser light, in which one of 
the mirrors acts as a nonlinear mechanical oscillator (see 
Fig. [T|) whose displacement is denoted by x. In addi- 
tion, the cavity medium is considered to be lossless, e.g., 
vacuum, and all optical losses (such as absorption and 
diffraction losses) occur at the mirrors. 

We refer the reader to the extensive body of litera- 
ture which exists for an in depth treatment of optical 
resonance cavities (see, for example, Refs. [42-46] and 
references therein). Here, we state the results which are 
needed in order to describe a simple optomechanical sys- 
tem. 

If the energy stored in the optical cavity in steady state 
reaches a local maximum at a; = xq, the intra-cavity 
optical power incident on a mirror can be written as |42l | 
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where F is the full width at half maximum parameter, L 
is the distance between two successive resonant positions 
of the micromechanical mirror, and 



where /pump is the power of the monochromatic light in- 
cident on the cavity, and Cre is the ratio of the resonant 
enhancement of the intra-cavity power. Note that for an 
empty cavity with metallic mirrors. 
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where A is the optical wavelength. In addition, the finesse 
of the optical cavity can be expressed as 

If the maximum mechanical displacement max|a;| is 
significantly smaller than F, a quadratic approximation 
for I{x) can be employed. In this case. 
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FIG. 1. (Color online) A general optomechanical resonance 
cavity. The left mirror is static. The right mirror is a me- 
chanical oscillator which can move in the direction parallel 
to the cavity axis [x direction). The cavity is pumped by a 
constant monochromatic light beam with the power /pump. 
The optical power circulating inside the cavity I depends on 
the actual position of the mechanical mirror, i.e., / = I(x). 
When the mechanical mirror is at rest, and no light is present, 
the mirror's position is denoted as a; = 0. The position of the 
mirror at which the optical power inside the cavity is maximal 
is called the spatial detuning and is denoted as xq. 



and 



I{x) w Jo + I'^x + -I'^x^ 
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where a prime denotes differentiation with respect to the 
displacement x. 

In this work, the optical power is assumed to follow 
the displacement without any lag. Namely, the optical 
response time is assumed to be much shorter than any 
other timescale in the system, including thermal relax- 
ation time and mechanical vibration period. 

The function I{x) in Eq. ([T]) can be represented by 
spatial Fourier series. 
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where 
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In practice, if the optical power I{x) changes relatively 
slowly with displacement, the series in Eq. ^ converges 
quickly (an example is shown in Fig. [T2|) . and Eq. ([4]) can 
be approximated as 



(6) 
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where fcmax }t J'- The exact expressions for Ck are derived 
in Appendix [Al 



B. Equations of motion 

We model the dynamics of the micromechanical mirror 
in the optical cavity by approximating it by a nonlinear 
mechanical oscillator with a single degree of freedom x 
operating near its primary resonance [43|. The mechani- 
cal oscillator's equation of motion is given by 

'^O .2 S 2 • 

X + —x + uj^x + asx + 73X X 

= 2/„cos(tJo + cro)t + F,p{x) + Fth(x), (7) 

where a dot denotes differentiation with respect to time 
t, X is the mirror displacement, wq is the original reso- 
nant frequency of the mirror, Q is the mechanical quality 
factor, ujjn is the momentary resonance frequency, whose 
dependence on ujq and other parameters will be discussed 
below, Q!3 is the nonlinear (cubic) elastic coefficient, and 
73 is the nonlinear dissipation coefficient. In addition, 
fm is the external excitation force, (jq is a small detun- 
ing of the external excitation frequency from ujq, F^p is a 
force resulting from radiation pressure, and -Fth is a force 
resulting directly from temperature changes in the mi- 
cromechanical mirror (such force can be attributed, for 
example, to thermal deformations [,23i] or buckling). 

Below, we consider external excitation frequency de- 
tuning (Tq to be small, i.e., ctq ^ ^^o- In addition, the 
mechanical quality factor is assumed to be large, i.e., 
Q> 1. 

It has been shown previously that nonlinear effects can 
play an important role in the dynamics of micromechan- 
ical systems [47l.l48|. In our case, we assume that the mi- 
cromechanical mirror behaves as a Duffing-like oscillator 
with positive nonlinear dissipation 73 > (i.e., the un- 
coupled autonomous mechanical system {fm = 0) is un- 
conditionally stable). Note that througliout this study, 
the mechanical nonlinearities are assumed to be weak, 
i.e., asx^ <C and 732;^ <^ oJo/Q- 

We assume linear dependence of the mechanical reso- 
nance frequency on the temperature: 



= ujo-/3iT-To), 



(8) 



Note that = cl^ because I{x) is real. 



where /3 is a proportionality coefficient, T is the effec- 
tive temperature of the mechanical oscillator, and Tq is 
the temperature of the environment. In the majority 
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of experimental situations, /3 is positive, i.e., heating of 
the micromechanical oscillator reduces its resonance fre- 
quency, while cooling increases it. 

In general, the nonlinear coefficients and 73 are 
functions of temperature similarly to Wm- However, due 
to the fact that the nonlinear terms are assumed to be 
small in Eq. ([7]) and the impact of their thermal variation 
is much smaller than that of , we regard the nonlinear 
mechanical coefficients as constants. The same is true for 
the linear dissipation coefficient ujo/Q. 

The time evolution of the effective temperature is gov- 
erned by the following equation: 



t = k{To^T)+7]I{x), 



where 



'had 
mCm 



(9) 



(10) 



the effective mass of the oscillator is denoted by m, h^^d 
is the radiation absorption factor of the mirror material. 
Cm is the mass-specific heat capacity of this material, rj 
is the heating rate due to interaction between the ma- 
terial of the mechanical oscillator and the light in the 
optical cavity, Tq is the temperature of the environment 
and K is the effective thermal conductance coefficient. In 
this simple approximation, the nonuniform temperature 
distribution due to localized radiative heating in the mi- 
cromechanical mirror is disregarded. 

In general, in addition to radiative heating term ril{x), 
Eq. ([9]) should account for heating due to mechanical 
damping. The heating power of this process can be esti- 
mated as 



Pc 



1 mojlA^ 
T 2Q ' 



where t = 2n/ujo is the time period of the mechanical 
vibrations, A is the amplitude of these vibrations, and 
all nonlinear effects have been neglected for simplicity. 
Comparing this heating power to the heating term in 
Eq. ([HI), we find that Pq is generally negligible if 



2Q 



< 1. 
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For example, for typical values of wq = 10^ sec" ^, 
m = 10"^^ kg, Q — 10^, and A — 1pm, we find that 
Pq « 8 X 10"^^ W. We compare this to the radiative 
heating by assuming that the radiation absorption factor 
/iiad of the micromechanical mirror is of order of several 
percents. It follows that if the optical power / in the 
cavity is approximately 10 nW or higher, the radiative 
heating is the dominant heating process. In practice, the 
optical powers that can have a significant impact on the 
system's dynamics and that are used in the experiments 
are of order of microwatts or higher, and, therefore, a 
term proportional to Pq is neglected in Eq. ^ . 



The formal solution of Eq. © is 



{kTq + ril{x)) e'^^d.T + T{t = 0) 



This can be shown to result in 



(12) 



(13) 







where the initial transient response term 
g-Kt ^rp^^ = 0) — To] has been dropped as insignifi- 
cant to the long timescale dynamics of the system. 

Using the fact that the energy and the momentum of a 
photon follow the relation f photon = cpphoton, where c is 
the velocity of light, we find that the radiation pressure 
force is 



Frp{x) = Vl{x), 



(14) 



where 



2 

mc 

and where light absorption by the micromechanical mir- 
ror has been neglected. 

Finally, we introduce a temperature dependent force, 
which acts directly on the micromechanical mirror. In 
practice, this thermal force can arise from several effects, 
such as a deflection of a bimorph mirror layer due to heat- 
ing, or a distortion due to internal stress [i^ [50l | caused 
by a non uniform heating of the mirror. The thermal 
force Fth is assumed to be linear in the temperature dif- 
ference T — Tq, i.e.. 



th 



9{T - To) =dr] I{x)e^^^-*UT, (15) 



where Eq. (|T3|) has been used. 

The equation of motion ([7]) can rewritten in a closed 
form as 

X + -I- [ojQ — p-qK^)]^ x + ai,x^ + 73X^i 

= 2/„ cos(ljo + ^n)t + vl{x) + eriK{I), (16) 
where we have defined the functional 



K{f) 



f{t)e''^''-*^dT. 



(17) 



Before application of the combined harmonic balance 
- averaging method to Eq. (jl6p . we conduct a stability 
analysis of the full dynamical system defined by Eqs. ([7]) 
and ^ in App. |B] There, it is shown that Hopf bifurca- 
tion is possible in the original system, and the necessary 
and sufficient conditions for this bifurcation are derived. 
These conditions will be shown below to be very similar 
to those found using the slow varying evolution equa- 
tions. 
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C. High thermal conduction limit 



where 



For the case where the characteristic therraal relax- 
ation time is much smaller than any other time scale 
in the system, namely (^qI^ and Q/ujq, the equation of 
motion ()16p can be significantly simplified. The mem- 
ory kernel e"^'^"*-' in Eq. ([T7]) can be replaced by a delta 
function S{t — t)/K, i.e., 

T-To = ^Iix). (18) 

K 

Consequently, the equation of motion (I16p becomes 



^ujot + (j), 



(21) 



.. Wo . 

X + —X + 



n 2 



X -\- a^x^ + j^x'^i 



Wo I{x) 

K 

ly-^] I{x) = 2f,ncos{ujo + ao)t, (19) 

K J 



It is easy to see that if the thermal relaxation rate in 
the system is fast compared to the mechanical resonance 
frequency, then the sole result of the coupling between the 
mechanical system and the optical cavity is the addition 
of nonlinear elastic terms proportional to I{x), I{x)x and 
I{x)'^x in the mechanical equation of motion (fT9]) . The 
mechanical dissipation terms proportional to ujo/Q and 
73X^ remain unchanged. 



and where Ai and (f> are the oscillator's amplitude and 
phase, respectively, and Aq is the static displacement. 
Here, it is assumed that the amplitude Ai and the phase 
(j) do not vary significantly on a time scale defined by 
ujq^ and, therefore, can be considered constant during 
a single period of the mechanical oscillation. This as- 
sumption is commonly referred to as the slow envelope 
approximation. 

The details of the averaging process used to derive the 
slow envelope evolution equations are given in App. [Cl 
Here, we state the main results. 

Assuming all the frequency corrections as well as the 
static displacement Aq to be small, we find that [see Ap- 
pendix [C] 



^0 



+ lasAl 



where 



and 



= ujQ - —Pq = ujq - Awo, 

K 



(22) 
(23) 



D. Finite amplitude oscillations analysis 

In general, in order for dissipative terms to occur in an 
equation of motion, some retardation in the displacement 
dependent force acting on the system is required fo', '23| . 
In our case, it is the memory kernel integral in K{I) in 
Eq. ([T6| that provides this retardation. In other words, 
the finite thermal relaxation rate k and the coupling of 
momentary mechanical resonance frequency a;„i to the 
optical power I{x) can be expected to result in changes 
in the effective linear and nonlinear dissipation of the 
micromechanical mirror [see Eqs. ([5]), ([HI), and (fT7|) [. 

It follows from the above discussion that a nontriv- 
ial dissipation behavior can be expected when the rate 
of thermal relaxation k is comparable to the mechanical 
resonance frequency We investigate the dynamics 

of mechanical oscillations with arbitrary amplitudes, i.e., 
oscillations with amplitudes that can be comparable with 
the wavelength of the light. The behavior of the optical 
power / as periodic function of the displacement x has 
been described in Sec. Ill AI 

In order to solve the equation of motion (jl6p . we 
make use of a combined harmonic balance - averaging 
method 

It can be expected that if all the nonlinear and optic 
related terms in Eq. (|16p are relatively small then the 
motion of the mirror is very similar to the motion of a 
simple harmonic oscillator, i.e., 



k=-k„ 



P„(Ao,Ai)- rcke^^^^^Jnh^k^), (24) 



where J„(z) is the Bessel function of order n. The term 
Awo represents a small mechanical frequency correction 
due to the averaged heating of the micromechanical mir- 
ror vibrating with an amplitude Ai. 

The evolution equations are (see App. ICj) : 



^Al - PiTi^^ [ 2/3 Ao + — 
4- \ Wo 
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Wo 



and 



^10 = - ( (To + Awo - ^Al + P2/3?] 



sine/), (25a) 



Ai 



Swo 



2wo"'' ■ ' '"''k^ +4w2 

loq) Wo 



Wo 



■ cos 0, (25b) 



where a new slow varying phase variable has been defined 
as (recall that the detuning cto is assumed small) 



x{t) ~ Ai^^ A\ cosip, 



(20) 



(Tot. 
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Equations together with Eqs. and con- 
stitute a coupled set of first order differential equations 
describing the time evolution of the slow envelope of the 
solution of Eq. ([TC]) . Now, we proceed to explore two 
important special cases of the system's behavior - the 
dynamics at small oscillation amplitudes and the steady 
state solutions corresponding to various limit cycles. 



E. Small amplitude oscillations limit 



The equations and (pS)) can be significantly sim- 
plified for small oscillation amplitudes and static deflec- 
tions, i.e., for Aq^i <C F. To this end, we denote the 
oscillation amplitude Ai as Aig in this section and sim- 
plify Eqs. (Uni) and ^ to [see also Eqs. ©]: 



AujQs = — /o, 

K 

lo 



Aos 



07] 



(26a) 
(26b) 

(26c) 



where the oscillation frequency fls and the static deflec- 
tion Aqs are independent of the oscillation amplitude Ais- 
In this limit, P„ can be represented by the lowest order 
terms in its Taylor series expansion, i.e.. 



Pn{Ai,) 



E 



and 

Aq w Aqs 



ujQK i / ^Ori\ ' 



(30b) 



Finally, Eqs. (|25)) can be simplified to 



Ms = -7^is - ^ sin( 

4 ujo 



Ais 



(cto + Alos) Au + 7 - — cos ( 
4 Wo 



where 



^ , 73 ^2 
2Q 2 



LOq 



2a;o 



(31a) 
(31b) 

li (32a) 



AuJs = Aujos - -z — ^Os 
ZUJo 



2wo -I- 



2^0 



and 



_ 3a3 /3ry 3K2-t-8a;g ^„ 
^ 2wo 2k k2 +4^2 . 

r-ll + 3r, J" 



(32b) 

(33a) 
(33b) 
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Using the fact that 
d"/(a;) 



E 



k=-k„ 



n+ \ \ L 



J—) c.. 



(27) 
(28) 



we can make the following substitutions for Po.1,2: 



Pn 



E 



k=~k„ 



■nk 



Ai 



lo 



- Til a2 



^max 7 i 



k=-k„ 



P2 



k=-k„ 



(29a) 
(29b) 

(29c) 



Consequently, the equations for H. and ^0 [Eqs. ([^5]) 
and ((22t , respectively] can be expanded up to the second 
order in Ais and first order in Aujqs resulting in 



Bn 
4k 



(30a) 



It is customary to rewrite the evolution equations (j3ip 
in a complex form by defining the complex amplitude 



Aise't 



Ais+jAiscj)] e 



„]<i> 



A\s cos t/i = flge 



j'("o+o'o)t , 



(34a) 

(34b) 
(34c) 



where c. c. denotes a complex conjugate. Using these 
definitions, the complex evolution equation reads: 

jas + (i7 - ctq - AijJs)as + (g + ir)a^a* = -j^, (35) 

Evidently, the coupling of a micromechanical mirror 
to an optical cavity introduces two types of terms into 
the complex evolution equation (1351) - linear terms pro- 
portional to /q and nonlinear terms of the third order 
proportional to /q . In addition, the autonomous part of 
the complex slowly varying evolution equation (/m = 0) 
consists of an approximated Hopf normal form of the 
original system |38l . [5l| , and is expected to yield condi- 
tions for self-excited limit cycles following either a sub 
or supercritical bifurcation determined by the sign of the 
cubic damping coefficient r. 
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III. SMALL OSCILLATIONS BEHAVIOR 



A. Linear and nonlinear effects in the dynamics of 
tiie small oscillations 



The hnear terms governing the dynamics of the 
micromechanical mirror considered here are given in 
Eqs. (|32[l. The parameter AuJs describes a smah addi- 
tional resonance frequency correction which arises from 
changes in heating and elastic nonlinearity due to small 
static displacement Aqs- In general, this correction can 
be considered small, i.e., Aw^ <^ ujq. In contrast, the 
linear dissipation coefficient 7 can undergo significant 
changes as function of the optical power, resulting in 
qualitative changes in the system's dynamics. 

An optical power dependent effective quality factor 
QeB can be defined by 



1 


27 


Qcff 




1 






Q 



277- 



-OJn 



(36) 



Note that from the experimental point of view, the def- 
inition of an effective quality factor given above is con- 
venient because QeS can be extracted directly from the 
small amplitude free ring down measurements of the mi- 
cromechanical mirror. In addition, Qcff is a function of 
/q. It follows that the local properties of I{x) in the vicin- 
ity of X = have a profound impact on the effective linear 
dissipation of the system. If the micromechanical mirror 
is positioned at the negative slope of the optical response 
curve, i.e., if /q < 0, and optical power is large enough, 
then the effective linear dissipation can be significantly 
reduced, resulting in extremely large ring down times, 
or even become negative. Alternatively, if the mirror is 
positioned at the positive slope, i.e. if /q > 0, a signif- 
icant increase in the effective dissipation, also known as 
"mechanical mode cooling", can be achieved (see the dis- 
cussion and references given in the Introduction section 
of this article) . 

The possibility of a negative linear damping suggests 
that the micromechanical mirror can develop self-excited 
oscillations. This mode of operation will be further in- 
vestigated in following sections. Here, we calculate the 
threshold conditions for the linear damping 7 to become 
negative, namely, the value /th of /max and the value XQth 
of Xq at the threshold. 

Neglecting all nonlinear terms and terms proportional 
to Aluq/ujo, the self oscillation threshold condition at an 
arbitrary value of xq, as can be derived from Eqs. (j32ap 
and (HI]), is 



^0 

2Q k2 



V 



—Io+"t;]Io = 0. (37) 



It should be emphasized that under the assumptions de- 
scribed above, this condition coincides with the exact 



Hopf criterion in Eq. (|B4[) found in App. |B]for the origi- 
nal dynamical system defined by Eqs. (O and ©. 

For a system in which the thermal force is dominant, 
the term proportional to h' in Eq. (|37p can be neglected. 
In contrast, if the radiation pressure impact is much 
larger than any heating induced mechanical forces, the 
term proportional to 6 can be neglected. By demanding 
that the threshold optical power is minimal, we find that 



a;oth 



r 

2\/3 

r 
'2/5 



< -^lyln 



(38a) 



and 



/th 



uj„ K +UI0 4r 
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This threshold is shown in Fig. [51 and in Fig. [3] together 
with different stability regions. 

In order to better illustrate the changes in the linear 
damping coefficient 7 due to coupling to an optical res- 
onance cavity, we choose a set of realistic parameters, 
which are given in Table HI and draw the resulting 7 co- 
efficient for a range of xq and /max values. The result is 
presented in Fig. [5] 

The coupling of the micromechanical oscillator to an 
optical resonance cavity does not only introduce linear 
contributions to the equation of motion, but has an im- 
pact on the nonlinear behavior of the system as well. The 
evolution equation (|35p is characteristic for a Duffing- 
type oscillator with nonlinear damping (s^ . |48| . The 
nonlinear coefficients in Eq. pSI) . i.e., q and r, are func- 
tions of the second derivative of the optical power I(x) 
with respect to displacement. 

It follows from Eqs. (|33p that if I{x) is convex near 
X — 0, namely Iq > 0, then the nonlinear elastic param- 
eter q is reduced (softening behavior), and the nonlinear 
dissipation is increased [see Eq. (j33bp ] if compared to the 
purely mechanical value r ~ 73/2. In contrast, if I{x) 
is concave in the vicinity of x = 0, namely Iq < 0, then 
the nonlinear elastic parameter is increased (hardening 
behavior), and the nonlinear dissipation is reduced. At 
optical powers high enough, the nonlinear dissipation can 
become negative, suggesting the existence of a large am- 
plitude limit cycle in the system (see Fig. 

Using Eqs. it can be can be shown that the effec- 
tive nonlinear corrections to the mechanical equation of 
motion discussed above change sign when 
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Xq ~ ± 
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4r 
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(39a) 
(39b) 
(39c) 



In this = is one of the inflection points of I{x). 
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parameter 



LJO_ 

Q 

73 



V 



L 

r 

To 



value 



20 X 10-^2 

160 
2.5 X 10^ 
3 X 10^* 
9 X 10^^ 
7.3 X 10^ 
O.OOlwo ~ 10^ 
7.5 X 10« 
325 
4.7 
0.775 
0.12L = 0.093 
77 
25 



units 



27 / 



kg 
kHz 



sec 

1 

sec 
rad 
sec K 

K 
sec W 
sec 
kg m 

N 
k^K 

pm 



K 



TABLE I. Values of parameters in a numerical example used 
to illustrate the results of Sections III Dl and III El All the 
values are of the same order of magnitude as those found in 
our experiments, which are reported elsewhere (4ll. |47|. 



The magnitude of nonlinear effects in this system 
strongly depends on the ratio between the thermal re- 
laxation rate k and the mechanical resonance frequency 
ujQ. At very fast thermal relaxation rates, the elastic coef- 
ficient is g — > 3a3/2ajo, and also r — >■ 73/2. As expected, 
the heating dependent nonlinear terms become negligible 
when loq, which is a special case of the general result 
discussed in Sec. Ill CI 

At low values of k, i.e., when the thermal relaxation 
time is significantly smaller than the mechanical oscilla- 
tion period 27r/r2s, care should be taken when applying 
the results of the previous section, because the require- 
ment that Aloqs ^ 1^0 can be easily violated, making 
the Eq. (I5T|) and all the results following it in Sec. Ill El 
inapplicable. 



B. Transient behavior 

In order to demonstrate the complex dissipative behav- 
ior of our system, we consider the non excited (/„ = 0, 




- -3 



X 10 

-4 



xo /T 

FIG. 2. (Color online) Linear dissipation coefficient 7 vs. the 
spatial cavity detuning xo and the optical power /max- For the 
positive values of 7, the non dimensional parameter drawn, 
2'y/Q.s, is equal to the reciprocal of the effective quality factor 
Qeff [see Eq. (|36p ]. In the area above the thick black line, the 
linear damping is negative (7 < 0), i.e., the solution a; = is 
no longer stable. Note that for positive values of the spatial 
detuning xq, the linear damping can greatly exceed the pure 
mechanical value, 1/Q = 0.4 x 10"^. 



CTq = 0) solution of Eqs. (j3ip . which can be written as 



Aujs - 



Ai\=0. 



(40a) 
(40b) 



Equation (|40ap is a regular Bernoulli differential equa- 
tion, which can be brought to a linear form by a standard 
transformation y = A^^ . The solution is 



AUt) = 



Ai,(0)2e-27t 



l+4^Ai,(0)2(l-e-27t)' 



(41) 



where the initial condition is Ais{t = 0) = Ais{0). Equa- 
tion (I40bp defines a small correction to the free oscillation 
frequency. 

Several interesting cases can be distinguished in 
Eq. (|4ip . Figure [3] summarizes all possible cases of linear 
and nonlinear dissipation as function of the initial dis- 
placement xq and maximal optical power in the cavity 

If the nonlinear dissipation coefficient r is positive, 
only finite stable solutions of Eq. (|4ip exist. If the lin- 
ear dissipation coefficient 7 is also positive, then the sys- 
tem decays almost exponentially to a single steady fixed 
point Ais = 0. The rate of decay at times t > 7"^ 
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FIG. 3. (Color online) Linear and nonlinear dissipation co- 
efficients in an optomechanical resonator and self oscillation 
thresholds. The black dotted vertical lines limit the area in 
which Iq < [see Eqs. (|39p ]. The solid blue line denotes a self 
oscillation threshold above which the effective linear damping 
is negative, i.e., 7 < 0. The dashed red line denotes the re- 
gion in which the nonlinear dissipation is negative (r < 0), 
and, therefore, the small amplitude limit cycle Alc given in 
Eq. (|42p is unstable, suggesting that the existence of an ad- 
ditional large amplitude stable limit cycle is possible. 



is approximately equal to the linear rate 27. This de- 
cay rate of the optomechanical oscillations can be either 
larger or smaller than the pure mechanical dissipation 
rate, wo/Q + 73^93, depending on the sign of /q (see also 
Fig.©. 

In contrast, if r > but 7 < then the system decays 
not to a trivial zero solution but to a stable limit cycle, 
whose radius in the plane of the complex slow changing 
amplitude [38] is given by 



\a\lc 



1 



.1 

r 



(42) 



The convergence to the limit cycle is again exponential. 
The result in Eq. (|^^ is correct only if Alc is sufficiently 
small, i.e., if the assumption Alc ^ T holds. The oscil- 
lation frequency of this limit cycle can be found from 
Eq. (j40bll . resulting in the following expression for the 
phase variable ^: 



(43a) 



The limit cycle frequency ujq — Aujg is similar to the 
one extracted from the local stability analysis of the full 
dynamical system given in App. [B]in the limit Alc ~> 
[see Eq. (|B5|) ]. 

Unlike the unconditionally stable cases described 
above, the result given in Eq. (|4ip can diverge in finite 
time if the nonlinear dissipation is negative, i.e., r < 0. 
The divergence occurs if the denominator in Eq. (|4ip be- 
comes zero. Here, two cases should be distinguished. If 



the linear dissipation is positive, i.e., 7 > 0, then the sys- 
tem will diverge only if the starting point Ais(O) > Alc- 
In other words, the limit cycle described in Eq. (j42|) ex- 
ists, but is unstable. If, however, both linear and non- 
linear dissipation terms are negative - the solution of 
Eqs. (|40ap unconditionally diverges. The general large 
amplitude analysis which is applicable in the last two 
cases has been presented in Sec. Ill Dl 

At this point, it is possible to give an estimate of the 
divergence time too by requiring that the denominator on 
the right hand side of Eq. (|4T|) vanishes, i.e.. 



resulting in 



(44) 



1 / 1 47 



The approximate divergence times according to Eq. (|44|) 
are shown in Fig. Note that when the absolute value 
of 7 is very low, the divergence time too can be very long 
if the starting point Ais{0) is close to the unstable limit 
cycle (for 7 > 0) or the origin (for 7 < 0). This behav- 
ior can be especially important if the system dynamics 
is simulated numerically, in which case extremely long 

tr?infiient times are imHpsirj^ble 
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FIG. 4. (Color online) Approximate divergence time too as a 
function of the initial amplitude Ais{0) [see Eq. (|44|) |. The 
variables are chosen so that the ajces are dimensionless (l7|too 
vs. VAi40)2|r/47|). Two cases are shown: solid blue line 
represents the case of positive linear dissipation (7 > 0), 
dashed black line represents the case of negative linear dis- 
sipation (7 < 0). The nonlinear dissipation is negative in 
both cases (r < 0). 
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IV. SELF-EXCITED OSCILLATIONS 

It follows from the stability analysis in the previ- 
ous Section and in App. [B] that a system governed by 
Eqs. spontaneously develops self-excited oscillations 
if 7 < 0, and can also start self-oscillating if driven far 
enough from the stable region near the origin in case 
7 > and r < 0. Here, we derive the steady state so- 
lutions of Eqs. (|25p in order to give semi-analytical esti- 
mations of the amplitudes of the steady limit cycles that 
exist in the system and their frequencies. 

For convenience, we rewrite Eqs. (|25ap . (|23|) and ((22|) 
for a steady state solution (i.e., Ai = 0) without external 
excitation terms below: 



^Lc/r 



n = ujo Pq = ujo- Awo, 

K 



(45a) 



Ao = 



and 



(45b) 



!^ + llAl + 2P,(37^ — 



Wo 



Ai 



LOq 



( 2/3 Ao -f — 0. (45c) 
^0 J 



The small frequency correction at a given steady state 
amplitude Ai can be found from Eq. (I25b[) . resulting in 



Pi 

Ai 



2/3 




which corresponds to Aujg at small amplitudes [see 
Eq. (M ). 

In order to illustrate the various possible limit cycles 
that can occur in a system whose parameters are given 
in Table m we plot the non zero solutions of Eq. (|45c|) for 
a representative range of the mechanical cavity detuning 
Xfj and the optical power /max in Fig. [5l 

As can be seen in Fig. [SJ a limit cycle with non zero 
amplitude always exists when 7 < 0, but only exists for 
the higher values of optical power when r < 0. This can 
be explained by the fact that when the nonlinear dissi- 
pation coefficient r is already negative but close to zero, 
the limit cycle amplitude given by Eq. (|42|) is extremely 
large, and the small amplitude analysis is inapplicable, 
as explained in Sec. IIIIBI In other words, Eq. (|45c|) can 
have only the trivial zero solution even when the nonlin- 
ear dissipation is negative, but still close to zero. 

The steady state solution of Eq. (|45cp for xq = — 0.5r 
is shown in Fig. [5] The zero solution is stable as long as 
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FIG. 5. (Color online) The stable limit cycle amplitude vs. 
the spatial cavity detuning xo and the optical power /max. 
Plotted is the non zero solution of Eq. (|45c[l . normalized by 
the width of the optical power peak, F. The parameters of 
the optomechanical system are given in Table (T] Similarly to 
Fig. |3] the linear dissipation is negative (7 < 0) above the 
solid blue line, and the nonlinear dissipation is negative (r < 
0) above the dashed red line. The thin dash-dotted magenta 
lines represent the three values of a;o at which the steady 
state amplitudes vs /max are plotted in Figs. [6l [T] and [8] 
The corresponding spatial cavity detuning values are xq /V — 
-0.5, -0.02, and +0.02. 



the linear dissipation is positive, and a small stable limit 
cycle develops when 7 becomes negative. 

It is interesting to compare the case above, in which 
the nonlinear dissipation is positive when the zero solu- 
tion loses stability (see Fig. [6]) , with a case in which, as 
the optical power increases, the nonlinear dissipation be- 
comes negative before the linear dissipation does. Such 
for Xq = — 0.02r is presented in Fig. [71 As can be 
seen in this figure, two stable solutions and one unsta- 
ble solution coexist in a bistable region, whose limits are 
marked by vertical arrows. This results in an amplitude 
hysteresis when the optical power /max or the spatial de- 
tuning Xq are swept. 

The linear damping in the bistable region is positive, 
and, therefore, the zero solution remains stable. In addi- 
tion to the zero solution, another large amplitude stable 
solution exists, because the nonlinear damping coefficient 
r is negative. At small amplitudes, the amplitude of the 
separatrix, denoted by the dashed red line, corresponds 
to the solution of Eq. p2)) , which is marked by large red 
dots. At optical powers high enough, the linear damp- 
ing becomes negative, the separatrix amplitude reaches 
zero, and the only remaining stable solution is the large 
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FIG. 6. (Color online) Steady state amplitude as function of 
the optical p O WG r Jm ax at xo = — 0.5r. The optomechanical 
system's parameters are given in Table |T] The plot corre- 
sponds to a cross section of Fig. O which is defined there by 
the leftmost dash-dotted magenta line. The large black dots 
are the estimations of the stable limit cycle for small ampli- 
tudes as given by Eq. (|42|l . 



amplitude limit cycle. 

The third typical configuratioii of limit cycles in this 
system is presented in Fig. |51 where a case for xq — 
-|-0.02r is shown. Here, the linear damping is uncon- 
ditionally positive, therefore, the zero solution is always 
stable. In addition, when the nonlinear damping is neg- 
ative, another couple of limit cycles can exist with finite 
amplitudes, an unstable one, acting as a separatrix, and 
a stable one. 

In order to complete the picture of the different limit 
cycles which are possible in the optomechanical sys- 
tem under study, the slow envelope velocity, Ai [see 
Eq. (|25ap ]. is drawn in Fig. [9] as a function of the am- 
plitude Ai at the bistable region shown in Fig. [S] 

Several features of Fig. [5] and Eq. (|25ap should be em- 
phasized. First, the stable finite amplitude solution S is 
separated from the stable zero solution O by the unsta- 
ble solution U. Second, the pair of fixed points U and 
S appear in a saddle node bifurcation when the opti- 
cal power is increased (in the case shown in Fig. [HI this 
bifurcation has already happened). Third, the positive 
mechanical nonlinear damping, i.e., 73 > 0, is prevalent 
at large amplitudes, driving the slow envelope velocity 
Ai to large negative values, and, therefore, preventing 
the existence of any other limit cycles with larger ampli- 
tudes. If the nonlinear mechanical effects are negligible, 
the system can become mult ist a ble, with several coexist- 
ing large amplitude limit cycles [20 , [23| . 



FIG. 7. (Color online) Steady state amplitudes as functions of 
the optical power /max at xo — — 0.02r. The optomechanical 
system's parameters are given in Table IT] The plot corre- 
sponds to a cross section of Fig. (5] which is defined there by 
the middle dash-dotted magenta line. The solid lines corre- 
spond to a stable steady limit cycle (black) and a stable zero 
solution (blue). The dashed red line corresponds to the un- 
stable hxed point (i.e., separatrix). The system is bistable 
in a certain range of optical powers, whose limits are marked 
by thin vertical arrows. The large red dots are the estima- 
tions of the unstable limit cycle for small amplitudes as given 
by Eq. (I42|l . As expected, the red dots align well with the 
unstable solution of Eq. (|45cP at lower values of Aj^c- 



V. 



NUMERICAL VALIDATION AND THE 
LIMITS OF ACCURACY 



In order to validate the analytical expressions derived 
above in Eqs. (I?5|l . and (P^ . we compare them to 

the results of the direct numerical integration of Eqs. ([T]), 
^ , and ^ . The values of all parameters used in the nu- 
merical simulation are given in Table m The value of the 
optical power I{x) in numerical simulations is calculated 
exactly, i.e., fcmax — 00. The numerical integrations were 
done using the Matlab software. 

The numerical results for the stable limit cycle am- 
plitudes at Xo = +0.02r are shown in Fig. [TOl together 
with the semi-analytical (i.e., slow envelope approxima- 
tion) results already presented in Fig.|Sl The comparison 
yields good agreement. 

The slow envelope approximation gives an estimation 
of the oscillation frequencies associated with large limit 
cycles [see Eqs. (|45ap and and their small vibra- 

tion limit [see Eqs. ((43|) ]. In Fig. [Til the free oscillation 
frequencies extracted from the numerical integration re- 
sults are compared with the semi-analytical results given 
in Eqs. (H5a|) and ^ for xq = -h0.02r. 
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FIG. 8. (Color online) Steady state amplitudes as functions of 
the optical power /max at xo = +0.02r. The optomechanical 
system's parameters are given in Table |T] The plot corre- 
sponds to a cross section of Fig. O which is defined there by 
the rightmost dash-dotted magenta line. The solid lines corre- 
spond to a stable limit cycle (black) and a stable zero solution 
(blue). The dashed red line corresponds to the unstable limit 
cycle (i.e., separatrix). The system is bistable above a certain 
optical power. The large red dots are the estimations of the 
unstable limit cycle for small amplitudes as given by Eq. (|42|l . 
As expected, the red dots align well with the unstable solution 
of Eq. (|45cp at lower values of A\,c. 



The limit cycle oscillation frequencies calculated us- 
ing Eqs. (|45ap and have a reasonable accuracy only 
when Aw < O.Iwq. This is due to the fact that we have 
neglected terms proportional to powers higher than one 
of Awo in Eqs. and ^2f2\ . This assumption of small 
frequency shift becomes increasingly inaccurate at high 
optical powers, as can be seen in Fig. [TT] 

In general, the linear expression in Eq. ([S]) is valid for 
small frequency corrections and for small temperature 
changes only. The accurate relation between the mechan- 
ical frequency and the effective temperature is usually 
more complicated, and strongly depends on the specific 
mirror configuration. For example, if a uniform doubly 
clamped beam with high internal tension is used as a 
mirror, its fundamental mode frequency can be approx- 
imated by a frequency of a fundamental harmonic of a 
pure string (52j 
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FIG. 9. (Color online) The slow envelope velocity, Ai [see 
Eq. (|25ap ]. as a function of the amplitude A\ at the bistable 
region shown in Fig. [S] The spatial optical cavity detuning 
is xo = -1-0. 02r, the optical power is /max = 55 mW, and 
the optomechanical system's parameters are given in Table [D 
In order to plot a dimensionless function, the velocity A\ is 
normalized by the characteristic fast mechanical velocity uioV. 
The stable zero amplitude solution O is denoted by a blue dot. 
The unstable limit cycle U is denoted by a red cross, and the 
stable finite amplitude limit cycle S is denoted by a black dot. 
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FIG. 10. (Color online) Numerical validation of the slow en- 
velope approximation results in Sec. Ill Dl The limit-cycle 
amplitudes as given by the solutions of Eq. (|45c|) are com- 
pared to the results of a full numerical integration of Eq. ([7|. 
The optomechanical system's parameters and the notation 
are similar to those used in Fig. [8] In addition, the initial 
conditions for the numerical simulations (green crosses) are 
shown, connected by thin green arrows to the final numerical 
solutions (green circles). 
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FIG. 11. (Color online) Numerical validation of the slow en- 
velope approximation results for frequency shifts from the 
mechanical frequency ujo- Note that the frequency shift is 
defined as ujo — Alj, i.e., the oscillation frequency is reduced 
when Au is positive. The semi-analytical limit-cycle and free 
small oscillation frequency shifts as given by the solutions 
of Eqs. (|45ap and (|46p are compared to the frequency shifts 
extracted from the results of a full numerical integration of 
Eq. ([7]). The optomechanical system's parameters are simi- 
lar to those used in Fig. [8] and are given in Tabled! Dashed 
lines represent the frequency shift which is solely due to the 
averaged heating (Eq. (|45a|) for a large limit cycle (black seg- 
ment on the right) and small vibrations near the origin (blue 
almost diagonal line across the figure) . Thin solid lines of the 
same colors represent the more exact solution, which incor- 
porates both Eqs. l|45ap and H46[l. The results of numerical 
simulations are represented by green asterisks. 



and Wstring is the string's angular vibration frequency, 
m and L are the mass and the length of the string re- 
spectively, S{T) is the temperature dependent total ten- 
sion in the string, 5*0 is the tension at T = Tq. i? is the 
Young's modulus, and a is the thermal linear expansion 
coefficient. Here, it is assumed that both the difference 
between the relaxed beam length and its actual length, 
and the change in the spring's tension due to heating are 
small. In addition, the Young's modulus and the thermal 
expansion coefficient are assumed to be constant in the 
relevant range of temperatures. One should also remem- 
ber that the notion of a single effective temperature T 
may not be sufficient to describe the thermally depen- 
dent mechanical behavior of a complex micromechanical 
structure. 

Another limit on the accuracy of the model described 
in Sec. lIIDl stems from the small nonlinearity assumption 
made in the slow envelope approximation [53]. Specifi- 
cally, only if the contribution of the nonlinear terms in 



Eq. (O is much smaller than the magnitude of the linear 
terms in the same equation, i.e., only if 

and 



then the harmonic solution assumption in Eq. (j20p to- 
gether with the averaging process used in Sec. Ill Dl are 
valid. 



VI. SUMMARY 

A coupling between an optical resonance cavity and a 
micromechanical resonator presents an interesting chal- 
lenge for building a simple yet comprehensive model, 
which is able to capture the complicated dynamics of the 
coupled system in a small set of relatively simple equa- 
tions of motion. In this work, we have created such a 
model for a low finesse optomechanical resonance cavity 
in which the elastic element is realized in the form of a 
vibrating nonlinear micromechanical mirror. 

The optomechanical cavity is assumed to be constantly 
pumped by monochromatic laser light. Due to the low 
finesse of the cavity, the optical response time is consid- 
ered to be very fast compared to the mechanical reso- 
nance frequency, and, therefore, the optical power inside 
the cavity can be described as an instantaneous function 
of the mirror's displacement [see Eq. ([T])]. Under these as- 
sumptions, we write a set of coupled differential equations 
which describe the mechanical and thermal dynamics of 
the system [see Eqs. ([T]) and (jH]), respectively]. 

The optical power influences the micromechanical mir- 
ror's dynamics both directly in the form of radiation pres- 
sure, and indirectly through heating. Radiative heating 
causes the mechanical resonance frequency to change ]see 
Eq. ([5])]. In addition, a direct thermal force can exist in 
a mirror in the form, for example, of a bimorph ther- 
mal actuation ]see Eq. (|15p ]. The important property of 
all heating dependent forces is the retardation that they 
introduce into the equations of motion, which results in 
changes in the effective dissipation in the micromechani- 
cal system. 

The micromechanical mirror itself is described as a 
Duffing-like weakly nonlinear oscillator with nonlinear 
(cubic) dissipation. The motion of the mirror can be 
approximated by a simple harmonic function with slow 
varying amplitude and phase. Averaging over a single 
"fast" period of mechanical oscillation results in a set of 
slow evolution equations for the slow varying amplitude 
and phase. These equations are given for the externally 
excited case in Sec. HlDl and for the case in which no ex- 
ternal excitation exists - in Sec. IIVI In addition, estima- 
tions of the oscillation frequency and the static deflection 
are derived in Sec. Ill Dl 
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Unfortunately, the full evolution equations for arbi- 
trary amplitudes do not have a simple analytical solu- 
tion. However, they do have a convenient semi-analytical 
closed form, and can be readily solved by any software 
designed for numerical calculations, such as the Matlab 
package used in this work. The solution of the first-order 
evolution equations requires significantly less computing 
power than the full numerical integration of the origi- 
nal equations of motion, which can be computationally 
prohibitive, especially in the case of low damping rates 
and very long transient times. One must bear in mind, 
however, that a slow varying envelope approximation of a 
general dynamic system may have a deficiency of missing 
additional nonlinear phenomena such as coexisting multi- 
stable limit cycles, quasi-periodic response (due to in- 
commensurate external and limit-cycle frequencies), ho- 
moclinic bifurcations and possible chaos. 

The evolution equations can be further simplified if 
the mechanical amplitude is small. It has been shown in 
Sec. Ill El that both linear and nonlinear terms originating 
from the optomechanical coupling can be found in the re- 
sulting small amplitude complex evolution equation psp . 
The changes in the effective linear and nonlinear dissipa- 
tion, which are functions both of the spatial cavity detun- 
ing and the pumping optical power, are most important 
[see Eqs. and (155)) ]. For example, if the spatial cavity 
detuning is negative, the effective linear dissipation can 
become negative at optical powers above a certain thresh- 
old, causing a small limit cycle (i.e., self oscillations) to 
appear. The threshold, the frequency, and the ampli- 
tude of these small self oscillations can be predicted with 
reasonable accuracy using the small amplitude approxi- 
mation [see Eq. (|42|) and Fig. [6]. These results coincide 
with the predictions of the stability analysis of the full 
dynamical system which is given in App. |B] 

Even when the linear effective damping remains posi- 
tive, a stable limit cycle with a large amplitude can co- 
exist with a stable zero solution in the region in which 
the nonlinear damping is negative. In such a case, a hys- 
teresis in the self oscillation amplitude is possible in the 
system when either the optical power or the spatial cavity 
detuning are swept back and forth. All the possible situ- 
ations leading to self oscillations have been summarized 
in Sec. |lVl 

Finally, we compare the results which are derived from 
the slow envelope evolution equations with the full nu- 
merical integration of the original equations of motion in 
Sec. El As expected, the semi-analytical results of this 
work are well-correlated with the full numerical integra- 
tion results as long as the major assumptions of the slow 
envelope approximation are satisfied. In other words, the 
validity of the majority of the results presented here de- 
pends on the assumption that all the optical dependent 
and nonlinear terms in the original equation of motion ([7]) 
are small. 

In our treatment, the dependence of the different terms 
in the equation of motion on the effective temperature of 
the vibrating mechanical element has the simplest, i.e.. 



linear, form. In general, the method of slow envelope and 
the averaging technique used in this study can be utilized 
in order to deal with more complex and more realistic re- 
lations between the heating and the oscillation frequency 
or the thermal force. In addition, further development of 
the ideas presented above may incorporate a treatment 
of large frequency changes due to heating and a depen- 
dence of additional parameters, such as nonlinear elastic 
coefficient and all mechanical dissipation coefficients, on 
temperature. 

Based on the theory presented here, an experimental 
study was conducted by us, which was reported else- 
where (4l| . A comparison between the experimental re- 
sults and the theoretical model developed in this article 
yields a good agreement. In particular, the quantitative 
theoretical model successfully predicted the experimen- 
tally measured changes in the linear effective damping, 
the cubic nonlinearities, the threshold of the self oscilla- 
tions, the frequency and the amplitude of the self oscilla- 
tions, and the resonance frequency of the micromechan- 
ical mirror under different conditions. The experimental 
study was done using micromechanical mirrors with two 
different geometries and material compositions. 
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Appendix A: Spatial Fourier series of a periodic 
optical power function 

In order to calculate an analytical expression for Ck 
in Sec. Ill Al we proceed as follows. We rewrite Eqs. ([T]) 
and (HI) as 



/(x) = 



hl„ 



1 + h — cos y 



where 



X - Xq 
L ' 



(Ala) 
(Alb) 

2 VL/ ' 

^ Cke^^^'^'^ ■ (Aid) 

Multiplying both sides of Eq. (jAla|) by 1 + h — cosy, using 
the fact that cosy = (e^'^ + 6^^^) /2, and separating terms 
corresponding to different harmonics, one finds 

(1 + h)l3k - + Pk+i) = I [J^'"^'^ : ^ J U . (A2) 
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FIG. 12. (Color online) Comparison between an exact peri- 
odic function describing the optical power in a resonator given 
in Eq. ^ (solid black line) and its truncated spatial Fourier 
series. The optical resonance width is F = O.OSL. Red dotted 
line corresponds to fcmax = 3. Blue dashed line corresponds 

to /i)max ~ 15. 



Note that Pk = P-k, and /3k are real because I{y) is a 
real even function. Assuming that P^ can be represented 
as 



(A3) 



Pk AnaxX*-^^ ^ 5 



where x and a are real, and substituting Eq. (jA3p into 
Eq. (jA2p for positive values of k results in 

- 2(1 + + 1 = 0. 

The solution which ensures series convergence by satisfy- 
ing the condition < a < 1 is 



a = l + h- y/{l + h)^ - 1. 



(A4a) 



The value of x can be found from Eq. (|A2[) for the case 
in which k = 0, giving 



Finally, Eq. (|Aldp gives 



(A4b) 



(A4c) 



It is straightforward to show that if the finesse is big- 
ger than unity, i.e., T is of order of ten or higher, the 
truncation error in Eq. © is negligible if fcmax ^ ^■ 

An example of several truncated Fourier series calcu- 
lated using Eqs. (jA4p for different values of fcmax is shown 
in Fig. [m 



Appendix B: Equilibrium analysis of the equations 
of motion 



In this section, we analyze the equilibrium position of 
the third order autonomous nonlinear dynamical system 
defined by Eqs. (O and ([9)) where the external exciting 
force is zero (/,„ = 0). 

By defining new variables p ~ x and AT = T — To, the 
equations of motion can be rewritten as 

X = p, (Bla) 

Wo 



732^ ] P ~ uj^x — asx +9AT + iyI{x), 



AT = -KAT + riI{x) 



(Bib) 
(Blc) 



where parameters defined in Sec. Ill Bi have been used. 

The equilibrium position of the dynamical system (i.e., 
the fixed point) is readily obtained by setting the veloc- 
ities (i.e., the left hand side of Eqs. (jBlf) ) to zero. This 
results in a transcendental function for the equilibrium 
displacement Aqs, 



n'iAos + aaAl - lyliAos) - 0ATo = 0, 
where the equilibrium temperature shift ATo is 

ATo = -HAos), 



(B2a) 



(B2b) 



and the equilibrium mechanical resonance frequency is 

^ujo- I3ATq. (B2c) 

In the limit of a very small equilibrium displacement 
Aqs ~ (i.e., the limit of very weak optomechani- 
cal forces), the Eqs. (jB2|) converge to the similar equa- 
tions dH]) derived in Sec. UTeI 

In general, multiple solutions of Eqs. (jB2[) may co- 
exist, corresponding to several stable and unstable fixed 
points under the same experimental conditions. How- 
ever, in the case in which the thermal frequency shift, 
the radiation pressure and the thermal force are all con- 
sidered small, the limiting case of Eqs. (pS)) predicts a 
single stable fixed point with a small static displacement 

^0. < r. 

Stability of the equilibrium is obtained via a local per- 
turbation of the system fixed point defined by Eqs. (jB2l) . 
resulting in a linear variation 







( 




= M 




Vat; 




^ AT 



where M is the Jacobian matrix of the first derivatives 
of the system functions given by the right hand parts 
of Eqs. (|B1[) . Thus, equilibrium stability can readily be 
obtained by evaluating the eigenvalues Ai, A2 and A3 of 
M, which satisfy: 

A^ -l-ciA^ + C2A + C3 = 0, 
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where 



ci 



(B3b) 

C3 = At {^ll + SasA^J - [nv + 77(2/3r!,Ao, + 9)] /'(Aqs), 

(B3c) 

and where a prime denotes differentiation with respect 
to the mechanical displacement x. 

Asymptotic stability of the equilibrium (i.e., Re{Ai} < 
0) is defined by positive coefficients and a positive sec- 
ond Hurwitz determinant, namely, > and A2 — 
(ciC2 — C3) > 0. Loss of equilibrium stability is defined by 
a zero eigenvalue (03 = 0) , or a Hopf bifurcation where 
the Jacobian matrix M has a pair of pure imaginary 
eigenvalues, i.e., Ai^2 = iiw^f. 

The zero eigenvalue condition C3 = can be rewritten 
in a differential form as 



(B3a) 



+ — I dI{A„s). 
K J 



This equation can be readily understood as a condition of 
equality between the thermally dependent nonlinear elas- 
tic force (left hand side terms) and the optomechanical 
forces (right hand side terms). This condition describes 
a saddle-node bifurcation, which can be reached for the 
case of larger optomechanical coupling than considered 
in this work. Note that the validity of the assumptions 
made in Sec. IIIBl especially the linear temperature de- 
pendence of the mechanical frequency and the thermal 
force, has to be carefully assessed in this case. 

The Hopf bifurcation, which implies that periodic limit 
cycle oscillations can occur near the bifurcation threshold 
can readily be shown to correspond to a zero second 
Hurwitz determinant, i.e., C1C2 — C3 = 0, with a positive 
Hopf frequency ujh = ^/c2- Using Eqs. (jBSp . we find the 
bifurcation threshold condition to be 



Wo ,73-2 



2Q ' 2'^^^^''^^ 

_ 2Q 2 ^0 



iyI'{Aos)~3a3Al 



I'iAos) 



Wq 

2Q 



-A, 



Os 



If we assume the mechanical dissipation, the nonlin- 
ear effects and the optomechanical coupling to be weak, 
namely, we assume the thermal frequency shift, the static 
displacement, the nonlinear and dissipation terms, the ra- 
diation pressure and the thermal force to be small, and, 
therefore, neglect all the small terms of the second order 



and higher, then the right hand side of Eq. (jB4|) van- 
ishes. In this limit, the Hopf bifurcation condition given 
in Eq. (jB4|) coincides with the condition 7 = discussed 
in Sec. inEl [see Eqs. ((5^ and 

Under the same assumptions, the Hopf frequency be- 
comes 



WH = \/c2 « Wo - Aws, 



(B5) 



where Auig is defined in Eqs. (|26ap and (j32bl) . This result 
coincides with the limit cycle frequency expression given 
in Eq. (|43p in the limit of vanishing limit cycle amplitude. 

We note that the Hopf bifurcation can either be super- 
critical or subcritical, culminating with stable or unsta- 
ble self-excited limit-cycle solutions which are discussed 
in Sec. HVl 



Appendix C: Averaging of the equations of motion 

Using Eqs. © and (PU)) . we write the optical power 
expression / as 



Cke 



j"27rf (Ao+Ai cosi/i) 



(CI) 



k=-k„ 



It is beneficial to use the Jacobi- Anger expansion 

00 

gjzcosj ^ Jo{z)+2j2j"Mz)cosn£_, (C2) 



where z and ^ are some real variables, and J„(z) is the 
Bessel function of n-th order. The optical power expres- 
sion given in Eq. (jCip can be rewritten as 



2 ^ Pn COS mp, 



(C3) 



where P„ are defined in Eq. ([M)) . 

Next, we proceed to write the integral in Eq. ((13]) ex- 
plicitly. Slow envelope approximation implies that the 
amplitude Ai, and the phase (f) do not undergo signifi- 
cant changes at timescales comparable to w^^. It follows 
that Ai and (p can be regarded as constants at timescales 



of order oJq and n , and terms involving K in Eq. ([16 



(B4) can be estimated using the approximate equality 



/(t)5(t - t)dT 



fit) / 9ir 
Jo 



t)dT, 



(C4) 



where g{T - t) is either e'^^^"*), e^>^±i^o){-r~t) 
g(K±j"2ajo)(T-t)^ /(t) is a function of slow varying terms 
Ai and <j), and all fast decaying terms in J g[T — t)dT 
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should be neglected. The result is 



K{cosnip) — / cosmpe'^^'^ '■'d? 
Jo 



1 



+ e 



K COS nip + nujQ sin nij) 



where f2 is defined in Eq. (j23p , and terms proportional to 
[firjKY have been neglected because the frequency cor- 
rection due to heating is considered small, i.e., l3rjK{I) <C 
LjQ. The term Awq can be identified as a small frequency 
correction due to the heating of the mirror averaged over 
one mechanical oscillation period. 

Equation (|C9I) can be further simplified by assuming 
the static displacement Aq to be small and using the weak 
nonlinearity assumption, i.e., a^A^ <C loq. giving rise to 
Eq. 1221). 

The remaining terms in Eq. (jCSP constitute the follow- 
ing relationship [see also Eqs. (|C3[) and (jCSp ] : 



^2 ^2|^2 



In order to solve Eq. (jl6p under the conditions de- 
scribed above, we use the harmonic balance method fol- 
lowed by the Krylov-Bogoliubov averaging technique , 
and require 

X = Aq + Ai cosip, 

X = —loqAi sinip, (C6a) 
X ~ —ujqAi cosTp — uJoAi simp — uJoAicf) cos ip . (C6b) 

It follows that 

AiCOSTp - Ai4>sintjj = 0, (C7) 
Introducing Eqs. (jC6|) into Eq. results in 

— ljqAi cosip — ujqAi sinip — ujoAicpcosip 
- sinip + [ljo - NK{I)f (Ao + Ai cos^A) 

+ [03(^0 + M cosip) - 73r2Ai smip] {Aq + Ai cosip)'^ 
= 2/„ cos(cjo + cra)t + vl + 9rjK{I). (C8) 



Ai sin ip + Ai(pcos ip — B, 



(CIO) 



where 



B 



— +73^^+4^2/3??- 



4 + \ cjo 



sini/j 



K 



803 



+ 2P2/3r/ 



+ 4a;o 



3a3 



2A?7^^ (2/3^0 + - )-2Pi- 

+ Cjg V Wo / Wo 



4wo 
cosV' 



Al 



cos(wo + (To)i + NST. (Cll) 



Here, NST denotes the non secular terms (i.e., higher 
harmonics) . 

The Eqs. (|C7[) and (jC10|) can be rearranged as follows: 



Ai — B sin-0, 
Ai(p ~ B cos-0. 



(C12a) 
(C12b) 



Averaging of Eqs. (jC12l) over one period of ip can be made 
under the assumption of slow varying envelope, namely: 



Collecting all non-harmonic terms in Eq. (jC8[) gives 
the expression for Aq: 



asAl +{n^ + ^asAlj Aq 

= ^P^Pv^^^+Po(-+'-^), (C9) 



27r 



Ai = — / Bsinipdip, 
27r Jo 



Ai^ 



27r Jo 



B cos ipdip. 



(C13a) 
(C13b) 



Substituting Eq. ([CTT|) into Eqs. ([CT3| yields the slow 
envelope evolution equations (|25p. 
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